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Abstract. We give a detailed description of the recently developed multi-domain spectral 
method for constructing highly accurate general-relativistic models of rapidly rotating stars. 
For both 'ordinary' and 'critical' configurations, it is exhibited by means of representative 
examples, how the accuracy improves as the order of the approximation increases. Apart 
from homogeneous fluid bodies, we also discuss models of polytropic and strange stars. 

Key words, stars: rotation — stars: neutron — gravitation — relativity — methods: nu- 
merical 

1. Introduction 

The structure and the gravitational field of relativistic, axisymmetric and stationary, uniformly 
rotating perfect fluid bodies is investigated in order to model extraordinarily compact astrophys- 
ical objects like neutron stars. The numerical calculation of these objects was the subject of 
papers by several authors, see Bonazzola & Schneider 1974, Wilson 1972 Butterworth & Ipser 
fT975lfT976l Friedman et al. fT986lfT989l Lattimer et al. fT990l Neugebauer & Herold fT992l Herold 
& Neugebauer [T9^2l Komatsu et al. U989al [T989b Eriguchi et al. [T9^4l Stergioulas & Friedman 
[T9951 Bonazzola et al. [T^9ll for reviews see Friedman fT^9"8l and Stergioulas 1998 

The idea to use a 'multi-domain spectral method' was introduced by Bonazzola et al. 
( 119931 [T998I I. In their 'BGSM-code', the space of physical coordinates is divided into several 
subregions, each one of them to be mapped onto a cross product of intervals. The physical field 
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2 Ansorg et al.: Highly accurate calculation of rotating neutron stars 

quantities are expressed in a spectral expansion with respect to all coordinates defined on the spe- 
cific intervals. If the interior region of the star is chosen to be exactly one of the domains, then 
it is possible to obtain representations of the field quantities that are smooth functions within the 
cross product of intervals. The spectral expansions then provide a very precise approximation of 
the field quantities. In particular, this choice for the domains circumvents the occurrence of the 
Gibbs phenomenon at the star's surface, which appears when non-smooth physical fields (such 
as the mass-energy density) are expressed in terms of a spectral expansion. 

The subject of this paper is a detailed description of our multi-domain spectral method (here- 
after 'AKM-method') which led to an accuracy of up to 12 digits for rapidly rotating, strongly 
relativistic homogeneous fluid bodies (Ansorg, Kleinwachter, & Meinel 2002 1. Its further devel- 
opment enabled us to calculate rotating toroidal bodies (the relativistic Dyson rings; see Ansorg, 
Kleinwachter, & Meinel 2003a I. 

The corresponding fundamental features of the AKM-method are as described above for the 
BGSM-code. However, the methods differ in the following: 

1. Instead of the three domains in the BGSM-code (two of them exterior to the star), we have 
only two domains, since we do not split up the region exterior to the star. 

2. As described in Bonazzola et al. ( 1998 1, the BGSM-code is an iterative scheme with each 
iteration step consisting of several procedures including the solving of nonlinear Poisson- 
like equations and the determination of an improved approximation of the star's surface. 
In the AKM-method, a large set of nonlinear algebraic equations for the unknown spectral 
coefficients corresponding to all field quantities and the unknown shape of the star's surface 
is simultaneously solved by a Newton-Raphson method. 

3. For the AKM-method, the restriction to only two domains, one of them exactly correspond- 
ing to the interior region of the star, can be maintained when the mass shedding limit is ap- 
proached. In this limit, the star is characterized by a cusp at its equatorial rim, thus requiring 
that the inner-most domain of the BGSM-code deviate slightly from the star's interior (see 
Bonazzola et al. [199 8, the resulting Gibbs phenomenon is limited since the displacement is 
small). 

We begin our description of the AKM-method with a review of the line element and Einstein's 
field equations, boundary, regularity and transition conditions as well as the resulting free bound- 
ary value problem. Following this, we provide an introduction to the method in question contain- 
ing a comprehensive description of the specific spectral field representations and the resulting 
set of nonlinear algebraic equations that ensures both the validity of the field equations within 
each domain and the transition conditions at the domains' boundary. In the subsequent section, 
we apply the AKM-method to 'ordinary' homogeneous, strange as well as polytropic stars (with 
a polytropic exponent r = 2) and display the improvement of the accuracy as the order of the 
approximation increases. Here, 'ordinary' stars means configurations that are not special in the 
sense that they neither rotate at the mass shedding limit, possess infinite central pressure, nor 
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show a considerable oblateness, but may be strongly relativistic as in the example given in table 
1 1 of Nozawa et al. fl998l and table 1 of Ansorg et al. d2002l 

Then we give representative examples concerning the excluded 'critical' configurations. At 
first, we consider homogeneous, strange as well as polytropic stars at the mass shedding limit. 
Following this, we study homogeneous configurations with an infinite central pressure. Finally, 
the last part of this section is devoted to highly flattened homogeneous bodies. 

In what follows, units are used in which the speed of light as well as Newton's constant of 
gravitation are equal to 1 . 

2. Line element, Einstein's field equations, and the free boundary value 
problem 

We use two different formulations of the line element describing the gravitational field of a uni- 
formly rotating perfect fluid body 1 . The corresponding Lewis-Papapetrou coordinates (p, ip, t) 
are uniquely defined by the requirement that the metric coefficients and their first derivatives be 
continuous at the surface of the body. 

Exterior to the star in question we write: 

ds 2 = e 2a (dp 2 + d?) + W 2 e- 2v (dip - cj dt) 2 - e 2v dt 2 (1) 

while for the interior we take metric functions valid in the comoving frame of reference: 

ds 2 = z w ' [e 2k '(dp 2 + df) + W 2 dp' 2 ] - e 2u '(dt + a' dip') 2 (2) 

Here, the only new coordinate is <p' — <p - Q.t where Q is the (constant) angular velocity of the 
star. 

Hence, we get the following transformation formulae: 
W-'e 2v ±(w-Q) = {We- 2U ' + a')~' , a = k' - U' (3) 
The exterior field equations following from the above line element read as follows 2 : 

a ( ;v = Iwvvj+<# (4) 

= 4(v p w p + v ( (jj ( ) (5) 
A 2 W = (6) 



1 The basic principles of the AKM-method rely on a rapidly converging Chebyshev representation of all 
physical and geometrical quantities within appropriate coordinates. The method is therefore applicable to 
arbitrary differentially rotating configurations (with some analytical rotation law). In this article we restrict 
ourselves to uniform rotation. In this case we may use the simpler field equations valid within the comoving 
frame of reference. 

2 Once v, ll> and W are known, the potential a can be determined by a line integral (with a — > as 

p 2 + -> co). 
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with the abbreviations 



a 2 = dl + dl 



,(./) 



a 2 + jw-HWpdp + w&y 



Since in the comoving frame the energy-momentum tensor reads 



T ik =(ji + p)iiu k + pg ik , u k = e~ u '6 k 4 



(7) 
(8) 

(9) 



where p is the total mass-energy density and p the pressure, the interior field equations assume a 
particularly simple form 3 : 



^ V) a' +A(a'U'+a'U') = 



A 2 W = 16 7re 



2(k'-U') 



Wp 



(10) 

(11) 

(12) 



As with the above potential a, the metric function k' can be determined via a line integration 
from the potentials U' , a' and W that follows from 

( kW PP - Wu) ) 



w p -w e 



w t w p 



w, 



Pi 



(13) 



' W[(U' p ) 2 - (U' ( ) 2 ] - \W- l z AU '[{a p ? - (a' ( ) 2 ] ^ 

2WU' p U' ( -\W- l & w a' p a' ( 
such that along the rotation axis p = the condition 

e = hm 

p^O p 

holds. Additionally, for a given equation of state, p = p(ji) or p = p{p), the relativistic Euler 
equations T' k k = yield p and p in terms of the metric potential U'\ 



(14) 



s u exp 



f 

Jo 



dp' 



e v ° = const. 



(15) 



o Kp')+p'_ 

Hence, for the exterior potentials v, a>, W as well as for the interior potentials U', a', W, particular 
systems of partial differential equations emerge. At the surface B of the star, the pressure p 
vanishes, which leads to a constant surface potential U' = Vo, see Eq. dl5l . If additionally the 
boundary B as well as the corresponding boundary values of the potentials a' and W were given, 
we would have to solve a particular interior and exterior boundary value problem of the respective 
field equations, completed by regularity conditions along the rotation axis (here, the p-derivatives 
of v, co, W/p, U' and a' vanish) and at infinity (here v — > 0, to — > and W - p — > 0). However, 
we have to deal with a free boundary value problem, where both the boundary B and the values 



See, for example, Kramer et al. ( 1980 1. 



Ansorg et al.: Highly accurate calculation of rotating neutron stars 5 

of a! and W at B are unknown, but have to be determined such that the normal derivatives of the 
potentials U' , a' and W behave continuously at B 4 . 

For a given equation of state, the corresponding solution depends on two parameters, e.g. the 
angular velocity Q and the gravitational mass M. Note that there might be multiple solutions 
corresponding to a specific prescribed parameter pair. For the description of the AKM-method, 
we consider at first the particular prescription of the pair (V"o, £2), but treat in a separate subsection 
the possible prescription of other parameter pairs. 

Together with the regularity conditions along the rotation axis, we assume that all metric 
potentials possess refiectional symmetry with respect to the equatorial plane ( — 0, leading to a 
vanishing f -derivative in this plane (see, for example, Meinel & Neugebauer l 1 9951 1. 

3. Description of the method 

A function, defined and analytic on a closed interval, can be represented by a rapidly converging 
Chebyshev-expansion. The spirit of the AKM-method is to use this property for all gravitational 
potentials, boundary values, and the unknown shape of the surface, which therefore need to be 
defined on appropriate cross products of intervals. The corresponding Chebyshev-coefficients are 
determined by a high-dimensional nonlinear set of algebraic equations that encompasses both 
field equations and transition conditions and is solved by a Newton-Raphson method. 

3.1. The mappings of the subdomains 

As already mentioned in the introduction, we divide the total space of physical coordinates into 
two subregions, an inner domain covering exactly the interior region of the star, and an outer 
one encompassing the exterior vacuum region. Both subregions are mapped onto the square 
I 2 = [0, 1] x [0, 1], which we realize by introducing a non-negative function y B defined on the 
interval / = [0, 1] that describes the surface of the body by 

B = {(p,0: p 2 = r 2 t, ? = r 2 v y B (t) , < t < 1} , 

(16) 

y B (0) = l, y B (l) = 0. 

Here r e and r p are the equatorial and polar coordinate radii of the body respectively. 
A particular example for the mapping in question is given by 

p 2 = r 2 e st, ? = r 2 p sy B (t), {s,t)el 2 (17) 

for the interior region and 

s s 

for the region exterior to the star. In this manner, the axis p = and the plane £ = are mapped to 
the coordinate boundaries t — and t — 1 respectively. Furthermore, the surface B of the body is 
4 It is a consequence of the field equations that k' is then also differentiable. 
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characterized by s = 1 . For the interior and exterior transformation, the point s — corresponds 
to the origin and to infinity respectively. 

Writing p 2 and ( 2 (and not p and £) in terms of the new variables s and t already takes the 
regularity condition along the rotation axis as well as the reflectional symmetry with respect to 
the equatorial plane into account. Indeed, for any potential that is analytic with respect to the 
variables s and f, it follows that the p-derivative at p — as well as the ^-derivative at £ = 
vanishes provided the above coordinate transformation is invertible there. The latter condition is 
only violated for s = 0. 5 

It turns out that the requirements of the regularity of the potentials (as functions of s and f) 
at t — and t = 1 replace a particular boundary condition here, that usually must be imposed. 
Similarly, the regularity of the interior potentials supersedes a boundary condition at the coordi- 
nate's origin. However, the asymptotic behaviour at infinity still must be considered, see section 

E21 

Note that for critical configurations we need to modify the above mapping, see sections l4.2.1l 
andlTOl 



3.2. The representations of the potentials and the surface 



For each of the gravitational potentials we use a specific Chebyshev-expansion that takes known 
boundary and transition conditions into account. In particular we know (r 2 = p 2 + ( 2 ): 



lim(rv) = -M , lim(r 3 w) = 2 J 



(19) 



lim r2(Wp _1 - 1) 



-2 



lim a'p 

p-»0 



limWp- 1 

p-»0 



< oo , 



< oo 



(20) 



where M and J are the gravitational mass and the angular momentum of the star, respectively, 
see also Eqs ( 1561 157> for an integral representation. Therefore we write outside the star: 

v = s[v B (t) + (s-l)H v (s,t)] (21) 



oj = s 3 [w B W + 0- l)HJs,t)] 

W (ext) = p (l + s 2 [W B (0 + (s - WwjaM]) 
and inside 



(22) 
(23) 



U' = V Q + (s- l)Hu>(s,t) (24) 
5 Note that the interior coordinate transformation introduced in section l4.2.3l is not invertible at the equa- 
torial rim of the star, but it is so at s = 0. 
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a' = p 2 [a' B (t) + (s-l)H ,(s,t)] (25) 

W (mt) = p [l + W B (t) + (s- l)H WM (s, 0] ■ (26) 

Here, the boundary values a' B , W B of the potentials a', W are expressed by the functions a^, Wb 
in the following manner: 

W B = p(W B + 1) (27) 

a B = P 2 a B . (28) 

The above one- and two-dimensional functions are expressed as limits of Chebyshev-expansions, 
e.g. 6 

v B (t) = lim vfHt) , (29) 



(m) 

V 



(0 = 2>^»>r w (2t-i), (30) 



H v (s,t) = lim H^\s,i) , (31) 



fl^Ci, 0=2 n^Tj-iV* ~ 1 ) r *-i( 2f - !)• (32) 

j,k=l 

Similarly, taking into account the representation of the boundary in d!6i . we write the boundary 
function y B as follows 7 

y B = (1 - [1 + rftg(.t)l (33) 



g{t) = lim g<- m Ht) , (34) 

m^co 



g {m \f) = Y J S ikM T k -x{2t-\). (35) 

k=l 

In the order m of our approximation scheme, we establish a nonlinear set of algebraic equations 
that determines the above coefficients of the m th Chebyshev-expansion. In the limit m — > oo, this 
set of algebraic equations is equivalent to the free boundary value problem in question, and the 
m th approximation becomes the solution. 



3.3. The nonlinear set of algebraic equations 

For a given equation of state, we specify the solution of our free boundary value problem by the 
prescription of a particular parameter pair. At first let us take ( Vo, £2); a more general choice will 
be discussed in section l3~5l 

6 The Chebyshev-polynomials are defined by Tj(x) = cos[/arccos(x)], x e [-1, 1]. 

7 In order to get m unknowns representing the surface of the star in the m th -order approximation discussed 

in section !331 we take the radii r c , r p and (m - 2) Chebyshev-coefficients for the function g. 
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We express the boundary values vb and wb in terms of (Vo, Q) and the functions a' B and Wb, 
see Eqs l|3}. This ensures the continuity conditions of the field potentials at the star's surface 8 . 
Hence, in the order m of our approximation scheme, we take the two-dimensional Chebyshev- 
coefficients 

Tj(jk;m) Tj(jk;m) u (jk\m) rj(jk;m) TT(jk;m) Tj(jk;m) 

Hy Maj , M Wext , ti v , , H a , , H WM (30) 

as well as the one-dimensional Chebyshev-coefficients 

(a' B f- m \W^ m \g {t ' m) (37) 
as independent variables. They build up a vector x (m) consisting of 

wJtotai = 6m 2 + 3m (38) 

components. The first 6m 2 components comprise all two-dimensional Chebyshev-coefficients 
while the following 3m-2 are the above one-dimensional Chebyshev-coefficients. The remaining 
two entries are filled by the values of r e and r p . 

We now describe in detail the components of a vector 

f(m) _ fOnVgCmH ^9) 

also consisting of m tota \ components that must vanish for the solution x (m) of the m -order ap- 
proximation. 

Given an arbitrary vector x (m) , we compute the Chebyshev coefficients corresponding to the 
first and second derivatives of the functions 9 

Tj(m) Tj(m) Tj(m) w (m) u (m) „(m) fAm 

and 

(a' B r\W B m \g^ (41) 

with respect to s and t. Together with the coordinate transformations (I17I18L we therewith find 
the first and second spatial derivatives of all gravitational potentials with respect to the coordi- 
nates p and f in our m th -order approximation, at an arbitrary grid point inside the domains (not at 
the origin or at infinity). So, we may fill the first 3m 2 entries of f (m) with the differences of right 
and left hand sides of the exterior equations (0}, evaluated at m 2 gridpoints (Sj, ft), j, k = 1 . . . m, 
corresponding to spatial points outside the star. Following the spirit of the spectral methods and 

8 The continuity conditions of the fields' derivatives will be part of the set of algebraic equations in 
question. 

9 Note that it is straightforward to calculate (i) the Chebyshev coefficients of a function from its values 
at the gridpoints (sj, t{) (or f* in the one-dimensional case), see Eq. 1421 . (ii) the value of a function at an 
arbitrary point inside or at the boundary of I 2 (or /) from its Chebyshev coefficients, and (iii) the Chebyshev 
coefficients of the derivative and the integral of a function from its Chebyshev coefficients, see Press et al. 
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in order to ensure a rapid convergence, we take for the sj and tj the roots of the m th Chebyshev 
polynomial, i. e. 

Sj = tj = cos 2 >0). (42) 

For the subsequent 3m 2 components of f (m) , we first compute the two-dimensional, m th -order 
Chebyshev coefficients corresponding to the interior function k'. This is done by determining the 
Chebyshev coefficients corresponding to the f-derivative of k' using Jl 3I > and again the coordinate 
transformation d!7i . and after that, by integrating with respect to the axis condition dl4> . see 
footnote [9] So, the 3m 2 entries in question can now be filled with the differences of right and 
left hand sides of the interior equations i ll Oi l, again evaluated at the m 2 gridpoints (Sj, ft) here 
corresponding to spatial points inside the star. 

The remaining 3m components of our vector f (m) are formed from the differences of the m th - 
order interior and exterior normal derivatives of the gravitational potentials v, u and W, evaluated 
at the m surface grid points (s = 1,^). Note that the coordinate transformations ( I17I18> are 
regular here, and thus the normal derivatives can easily be computed using the shape of the 
star that is incorporated in x ( '"K The interior normal derivatives of v and a> follow from the 
transformation formulae Q and the interior potentials. 

In this manner we get in the m th approximation order a nonlinear set of m tota i = 6m 2 + 3m 
algebraic equations 

f (m) (x (m) ) = (43) 
that is solved by a Newton-Raphson method, see section l3~4l 



3.4. The Newton-Raphson method and the initial solution 



In the Newton-Raphson method, the zero of a nonlinear set of algebraic equations of the form 
J43I is determined iteratively, 



Y (m) _ Y (»0 _ 
A n ~ A «-l 



Q f ('«) 




dx (m) 


JM 



f(m) ( (mK 

v n-V ' 



(44) 



requiring an initial which must already be sufficiently close to the final solution 



lim„^co xj, m \ The Jacobi matrix in the Eq. (1441 is determined approximately using (e <K 1) 



$f(»l) 




5x (m) 





AB 



I[fr(x^ + eej ,)-ff ( X «-ee B )l. 



(45) 



Here the subscripts A and B denote the corresponding element of the Jacobi matrix and the vector 
f (m) , and t B is the B th unit vector, (Cb)a = Sab- 

There are various possibilities for obtaining an initial solution. For example, one may start 
from the static solution characterized by Q. = 0. Here the corresponding field equations turn into 
ordinary differential equations with respect to the radial coordinate r, and these can be solved 
e.g. by using a Runge-Kutta-method. Taking this solution for the initial x^™', one may now grad- 
ually increase the parameter Q and thus eventually explore the whole parameter space. Another 
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possibility is to start with a Newtonian solution (e.g. a Maclaurin spheroid). Proceeding into the 
relativistic regime comes about by increasing the absolute value of Vq. 

In our treatment we favoured the latter initialization, for it already provides highly distorted 
bodies. Moreover, we calculated configurations with a particular equation of state by starting 
from a constant mass-energy density profile and continuously moving to the desired equation of 
state. 



3.5. Arbitrary parameter prescription 

With a slight modification of our nonlinear set of equations described in section l3~3l we are able to 
take various different parameter prescriptions into account. The idea is to add the quantities Q. and 
Vo to the vector x (m) , resulting in m tota i = 6m 2 + 3m + 2 unknowns from now on. Simultaneously, 
we add two equations to the nonlinear set representing exactly the desired parameter relation 
for the solution in question. This can be done since all physical quantities concerning the final 
solution are now contained in the vector x (m) . 

For example the potential U' c at the origin reads 10 

U' c = Vo-Hu,(0,0), (46) 

which is directly connected to the central pressure p c , see Eq. d!5l >. Likewise, M and J can be 
expressed 

M = -r p [v B (0)-# v (0,0)] , (47) 
J = ^[ WB (0)-/U0,0)]. (48) 

Also the prescription of a parameter ji is possible which controls the distance of a configuration 
to the mass shedding limit: 

dvu in the mass shedding limit 

£=— £ (t=l) = \ (49) 
I 1 for an ellipsoidal shape 

Similarly, one can prescribe more complicated expressions such as the baryonic mass Mq, which 
is defined by an integral over the interior field quantities. 

Any two conditions of this kind (of which the above ones are just examples) can be taken 
and added to the system of nonlinear equations. The corresponding parameters (here p c , M, J,/3 
or Mo) must then be prescribed. In this paper we concentrate on the pair (p c , r p /r e ) and only take 
iPcP) in order to place ourselve exactly on the mass shedding limit, see section l4.2.1l 



10 Note that the function H^\s = 0, /) tends to a constant in the limit m — » oo. Similar properties hold for 
all functions listed in PUl . see also section l331 
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3.6. Regularity and uniqueness 

As already depicted in section ETTl the AKM-method is characterized by the fact that some of the 
usual boundary conditions are replaced by regularity requirements. Moreover, if for the moment 
we only consider the functions 

H v , H^, Hw,ext,Hu> ,H a > ,Hw; m t (50) 

and treat the quantities r e , r p , £2, Vo as well as the functions a' B , Wb and g as if they were given, we 
obtain specific partial differential equations valid in I 2 , and particular boundary conditions at any 
edge of I 2 are not required for any of the functions listed in J50i . Nevertheless, the solution of this 
system of equations is uniquely determined if we require regularity with respect to all functions. 
A similar situation can be studied in the one-dimensional case, e. g. the equation ( ' = d/dt) 

f(l - tjh + 2(1 - 2t)h - 2h + 2 = o [f(l - t)h\" = -2 (51) 

possesses only the solution h = \ which is regular within /. 

The above approximation scheme sorts out the non-regular solutions since it is based on 
Chebyshev-expansions. It moreover ensures known, additional properties of the functions \5 01 at 
s — 0, e. g. Hw{s — 0, t) — const. Note that these properties are approached as m — > oo. 



4. Representative Examples 

4.1. Ordinary stars 

At first we apply the AKM-method to three models of homogeneous [Eqs (I52H . polytropic [with 
a polytropic exponent F = 2, Eqs d53i l. and strange stars [Eqs ( I54H . In particular, we prescribe 
the corresponding equation of state in the form fi = fi(p) and find from Eq. dl5> the relation to 
the interior potential U': 



fi (p) = fi = const 







(52) 



fl(p) = p+ yfp/K 



p = (e v °- u ' - if /(4K) 
p = (e 2 W.-f)_i)/(4Jr) 



(53) 



H(p) = 4B + 3p 



P = 



(54) 



H = B(l+3e 4 ^-^) 

Here, K and B are the polytropic constant and the MIT bag constant, respectively 11 . Note that 
for the application of the AKM-method (strictly speaking, only for its rapid convergence) it is 
11 For a description of the equations of state corresponding to polytropic and strange star matter, see 
Tooper 1 1965} and e.g. Gourgoulhon et al. 1 1999 1 respectively. 
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m=6 m=8 m=10 m=12 m=14 m=16 m=18 m=20 m=22 

Pc 1 

r p /r e 0.7 

Cl 1.41170848318 1.9e-04 1.3e-05 7.8e-07 2.9e-08 8.5e-10 4.6e-ll 3.0e-12 1.3e-13 8.0e-15 

M 0.135798178809 1.8e-04 3.5e-06 5.9e-08 3.4e-09 3.8e-10 2.6e-ll 8.5e-13 3.3e-14 6.8e-15 

R alc 0.345476187602 2.0e-08 1.5e-06 1.7e-08 1.8e-09 4.2e-ll 1.8e-ll 1.6e-12 l.le-13 1.3e-14 

/ 0.0140585992949 8.7e-04 6.8e-05 3.7e-06 1.2e-08 1.2e-08 6.8e-10 8.4e-12 3.5e-12 2.0e-13 

Z p 1.70735395213 3.2e-05 6.5e-06 2.4e-07 3.6e-09 4.6e-10 9.1e-12 7.1e-13 1.7e-13 1.6e-14 

GRV2 7.5e-05 3.9e-06 3.9e-07 2.2e-08 8.9e-10 4.2e-ll 3.1e-12 3.0e-13 7.7e-14 

GRV3 1.2e-05 7.5e-06 1.2e-07 2.9e-08 1.4e-09 3.5e-ll 1.3e-12 1.8e-13 6.5e-14 

\l-MJM oat \ 2.8e-04 4.9e-06 1.9e-07 l.le-08 4.1e-10 3.4e-12 1.5e-12 4.2e-13 2.3e-13 

|1 - Jin/^outl 1.2e-03 7.0e-05 4.1e-06 5.0e-08 l.le-08 7.1e-10 4.6e-12 2.9e-12 l.le-13 



Table 1. Results for a constant mass-energy density model (fi - hq) with p c = 1, r p /r e = 0.7. 

— 1 /2 - 1 /2 - 1/2 

1 1 ere. p c = p c /p.Q, Q = Q/yU , M — M/j. q , R c i rc = R c i rc iu and J = Jpu are normalized values 
of the physical quantities, see Eqs (1561 I57> . Apart from the virial identities GRV2 and GRV3 
in the m th order approximation, the columns 3-11 display the relative deviation of the specific 
quantity in the m th order approximation with respect to the numerical result obtained for m = 24. 
The quantities M- m , 7; n and M out , 7 out refer to the corresponding numerical values resulting from 
( 1561 157> and dl9> respectively. 







m=6 


m=8 


m=10 


m=12 


m=14 


m=16 


m=18 


m=20 


m=22 




l 






















0.834 




















Q. 


0.4004385709 


l.le-03 


9.0e-05 


7.3e-06 


6.4e-07 


6.4e-08 


7.2e-09 


8.6e-10 


l.le-10 


1.3e-ll 


M 


0.1605611357 


4.2e-04 


5.7e-06 


2.5e-06 


3.2e-07 


3.7e-08 


4.1e-09 


4.7e-10 


5.4e-ll 


5.9e-12 


^circ 


0.6794279802 


5.7e-04 


6.0e-05 


5.4e-06 


5.1e-07 


5.2e-08 


5.9e-09 


7.0e-10 


8.6e-ll 


9.8e-12 


./ 


0.009491087857 


9.4e-04 


3.3e-05 


8.5e-06 


l.le-06 


1.2e-07 


1.3e-08 


1.4e-09 


1.4e-10 


1.2e-ll 


z P 


0.4580590747 


1.7e-03 


8.7e-05 


4.6e-06 


2.6e-07 


1.9e-08 


1.8e-09 


2.0e-10 


2.6e-ll 


3.1e-12 


GRV2 




2.6e-04 


8.1e-06 


4.8e-07 


7.5e-08 


1.3e-08 


1.8e-09 


2.4e-10 


3.2e-ll 


4.3e-12 


GRV3 




5.5e-05 


2.9e-06 


2.7e-07 


1.5e-08 


1.2e-09 


2.0e-10 


3.9e-ll 


7.1e-12 


1.3e-12 


\l-M m /M out \ 




1.2e-04 


6.5e-06 


5.5e-07 


5.5e-08 


5.5e-09 


5.3e-10 


4.5e-ll 


2.4e-12 


3.3e-13 


1 1 — Jinf Jovtl 




2.7e-04 


3.7e-05 


4.0e-06 


4.2e-07 


4.1e-08 


3.6e-09 


2.5e-10 


5.6e-13 


4.9e-12 



Table 2. Results for a polytropic model (polytropic exponent T -2, polytropic constant K) with 
fi c = 1, r p /r e = 0.834. Here, fi c = Kfi c , Q = K l/2 Cl, M = R- l/z M, R circ = K- l/2 R circ and 
J = K l J are normalized values of the physical quantities, see Eqs ( 1561 \57i . For the meaning of 
the quantities listed in the columns 3-11, see Table 1. 
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m=6 m=8 m=10 m=12 m=14 m=16 m=18 m=20 m=22 

Pc 2 
r p /r e 0.5 

3.4304996 2.1e-04 1.5e-05 4.0e-07 1.6e-07 6.2e-08 1.9e-08 5.3e-09 1.5e-09 3.4e-10 

M 0.035510326 5.3e-03 5.9e-04 6.4e-05 7.2e-06 9.2e-07 1.4e-07 2.6e-08 5.8e-09 1.3e-09 

R alc 0.14117783 3.3e-05 3.2e-09 1.8e-07 8.4e-08 3.0e-08 9.4e-09 2.9e-09 8.2e-10 2.0e-10 

/ 0.0011024838 1.6e-03 2.8e-04 5.4e-05 8.5e-06 1.4e-06 2.8e-07 6.4e-08 1.6e-08 3.5e-09 

Z p 0.72634557 8.6e-05 1.8e-05 5.2e-06 1.3e-06 3.2e-07 8.4e-08 2.4e-08 6.3e-09 1.5e-09 

GRV2 1.2e-04 4.8e-06 4.1e-07 4.0e-08 3.9e-09 3.9e-10 3.8e-ll 3.2e-12 4.7e-14 

GRV3 6.7e-04 3.9e-05 2.8e-06 2.4e-07 2.1e-08 1.8e-09 1.5e-10 l.le-11 5.2e-14 

|l-M m /M 0Ut | 5.9e-03 6.0e-04 6.2e-05 6.5e-06 6.9e-07 7.4e-08 7.9e-09 8.4e-10 8.9e-ll 

|l-/in// ou tl 2.2e-03 2.7e-04 4.6e-05 6.0e-06 6.5e-07 8.1e-08 l.le-08 1.2e-09 1.5e-10 



Table 3. Results for a strange star model (MIT bag constant B) with p c = 2, r p /r e = 0.5. Here, 
p c = B l p c , D. = B 1/2 D., M = B 1/2 M, R circ = B l/2 R cilc and J = BJ are normalized values of 
the physical quantities, see Eqs ( 1561 157> . For the meaning of the quantities listed in the columns 
3-11, see Table 1. 



necessary to have analytic dependencies p = p(U') and p = p(U'), in particular at U' = Vo- For 
the equation of state, p. = p(p), this requires that 

^//Pl^/IM)] (55) 

where / is some function which is positive and analytic when its argument vanishes, and N is 
some non-negative integer. Apart from the homogeneous and the strange star model, this condi- 
tion is met for polytropic equations of state with a polytropic exponent F = 1 + 1 /TV when N is a 
non-negative integer (as in the case above where N = 1). In order to treat more general equations 
of state, one needs to consider several layers inside the star, with each one of them characterized 
by a specific equation of state. The outermost one of them again must meet the above require- 
ment. The consideration of several layers leads to a corresponding number of subregions into 
which the interior domain needs to be split. 

In Tables 1 to 3 one finds numerical values of several physical quantities, for a specified 
configuration with prescribed central pressure p c (equivalently, for non-homogeneous models, 
we may prescribe the central mass-energy density p c , see Table 2) and radius ratio r p /r e . The 
angular momentum J, gravitational mass M, equatorial circumferential radius R c i rc and the polar 
redshift Z p are given by : 

J = -In J (p + p)a'e 2k '- 2u 'Wdpd£ (56) 
M = 2Q7 + 2n j(p + 3p)s 2k '- 2u ' Wdpd( (57) 
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Rckc = e-MwVl -u 2 \ i 



0) 

(58) 



Z p = e- v ° - 1 



with u - -W~ l a'e 2u ' . Note that the above integrals extend over the space of p- and ^"-coordinates 
covering the interior of the body 12 . 

A first test of the accuracy of a solution determined numerically is the comparison of the 
calculations of M and J from the exterior fields [see Eqs dl9H with those from the above integral 
representations (1561 157> . A further check is given by the general-relativistic virial identies GRV2 
and GRV3, derived by Bonazzola & Gourgoulhon 1994 and Gourgoulhon & Bonazzola ll994l 
As a consequence of the field equations, they identically vanish for an exact analytic solution 
corresponding to a stationary and asymptotically flat spacetime. Particularly, for our rotating star 
models they read 

GRV3 = \l-h/h\, GRV2 = \l-Jr/J 2 \ (59) 
with : 



/, = 4tt 



3pVl-M 2 + (fl + P )- 



e 2k '- w 'Wdpd( (60) 



= Je" L [(Vv) 2 - ivaV/? - ^wV 4 " (Vco) 2 ~ (l - e 2 ^) (or, - i T/pJJ dpd£ (61) 



8tt \p + (p + p) 



u 2 



1 — u 



e 2k '- 2u 'dpd£ (62) 



J 2 = | |(Vv) 2 -^V 4v (V W ) 2 



dpd{ (63) 



(we use the abbreviations u = -W l a'e 2u , pe n = We v ). The Nabla-operator has its usual 
meaning, in terms of the coordinates p and f , i. e. for any two functions a and b 

(Va) 2 = a 2 + a 2 , VaVb = a p b p + a^b(, (64) 

The above integrals are taken over the whole space of p- and ^-coordinates (for I\ and J\ this 
reduces again to the interior region of the star since both pressure p and energy-mass density p 
vanish outside the body). 

Apart from the values for the above physical quantities with the accuracy that was reached in 
the 24 th approximation order, we provide in Tables 1-3 the improvement of the accuracy as the 
order m is increased. Also given are the corresponding numerical values of the general-relativistic 
virial identies GRV2 and GRV3 as well as those of the relative deviations concerning the integral 
and far-field representations of M and /. 

We note generally an exponential convergence of the numerical solution as the order m in- 
creases. This is a common feature of the spectral methods. However, the star's field quantities 
may vary in their 'smoothness' resulting in a variably rapid convergence. For example, the con- 
vergence of the numerical solution corresponding to the homogeneous star (Table 1) is much 
12 The quantities M and J can also be taken from the exterior fields v and a>, see Eqs ll9l 
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Fig. 1. Meridional cross-section of the homogeneous mass-shedding configuration specified in 
Table 4 (with the axes scaled identically). The dashed curves indicate the boundary of the corre- 
sponding ergo-region. 

faster than that corresponding to the strange one (Table 3). In a sense, the strange star is closer 
to the mass-shedding limit (here /3 ~ 0.39 while for the homogeneous model j3 » 0.84) and 
moreover more flattened. 

The models in Tables 1 and 2 have been calculated by Nozawa et al. ( 1998 ). Note that for the 
polytropic model there is a steeper gradient of the pressure as a function of the radial coordinate 
r, e.g. within the equatorial plane. In order to take this property into account, we used for this 
model the following slightly modified interior coordinate transformation 

p 2 = r e V(s;c s )f, ? = r*(ris;c s )y B (t), (s,t)el 2 (65) 
with 

1 -c s 

<t(s;c s ) = - s (66) 

1 - c s s 

and the constant parameter c s = 0.6 (for the other models we took c s = 0). Minor modifications of 
this kind, specially suited to the particular problem in question, may accelerate the convergence, 
see below for further examples. 



4.2. Critical stars 

4.2.1 . Stars at the mass shedding limit 

The endpoint of a sequence of rotating stars is often marked by a mass shedding limit. It is of 
particular interest since specific physical quantities such as the angular velocity reach maximal 
values there. A highly accurate determination of this limit is therefore desirable. 

The mass shedding limit is reached when the angular velocity Q of the star attains the angular 
velocity of test particles moving on a prograde circular orbit at the star's equatorial rim. For the 
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Fig. 2. First derivative of the function g = g(t) with respect to t for the homogeneous Newtonian 
mass-shedding configuration (A) of table 2 in Ansorg, Kleinwachter & Meinel ( 2003b I. The 
numerical methods explained in section 3 ibid, have been carried out up to the approximation 
order JV = 80. 







m=10 


m=16 


m=22 


Pc 


1 








r p /r e 


0.4918 


7.2e-04 


1.2e-04 


1.6e-05 


Q 


1.6588 


6.1e-05 


5.5e-06 


5.8e-07 


M 


0.1623 


7.2e-05 


3.8e-06 


6.1e-07 


"circ 


0.4041 


3.2e-05 


7.1e-06 


1.2e-06 


J 


0.02431 


1.8e-04 


9.4e-06 


1.6e-06 


z P 


2.6685 


1.2e-04 


1.6e-05 


2.1e-06 


GRV2 




1.2e-04 


1.6e-05 


4.3e-06 


GRV3 




1.3e-04 


1.8e-05 


4.9e-06 


\l-M-JM out \ 




2.8e-04 


4.9e-05 


1.4e-05 


U — J'ml >/oiitl 




4.2e-04 


8.1e-05 


2.3e-05 



Table 4. Results for a constant mass-energy density model rotating at the mass-shedding limit 
with p c — 1 . For the meaning of the quantities, see Table 1 . 

p-derivative of the field quantity U' it follows: 

f/;(r e ,0) = 0. (67) 

Moreover, a cusp at the surface occurs (see Fig. 1), which corresponds to a vanishing mass- 
shedding parameter/?, see Eq. ( 1491 . 

Numerical investigations of a homogeneous Newtonian configuration rotating at the mass- 
shedding limit suggest that the surface function g becomes singular in higher derivatives at this 
limit, see Fig. 2. This causes a similar singular behaviour of all gravitational potentials, and we 
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m=10 


m=16 


m=22 




0.34 










0.5845178 


9.1e-05 


9.7e-07 


3.6e-08 


Q 


0.3770150 


2.0e-04 


7.0e-07 


3.8e-09 


M 


0.1883522 


8.7e-05 


9.7e-08 


l.le-09 


"eirc 


1.0920220 


3.4e-05 


3.0e-07 


2.6e-08 


/ 


0.02023980 


3.0e-05 


1.4e-06 


9.5e-09 


z p 


0.4035809 


2.9e-04 


6.0e-07 


2.3e-09 


GtfK2 




1.5e-05 


8.7e-08 


1.6e-09 


Gi?V3 




5.5e-06 


8.4e-08 


1.3e-09 


ll-Mjn/AU 




6.7e-05 


2.7e-07 


3.2e-09 


|1 — «An/^outl 




2.8e-04 


1.6e-06 


1.8e-08 



Table 5. Results for a polytropic model (polytropic exponent F = 2) at the mass-shedding limit 
with fi c = 0.34. For the meaning of the quantities, see Table 2. 







m=10 


m=16 


m=22 


Pc 


3 








r p /r e 


0.4713 


6.0e-04 


1.4e-04 


1.9e-05 


Cl 


3.6505 


4.0e-05 


6.3e-06 


6.0e-07 


M 


0.03719 


1.5e-04 


9.4e-06 


1.3e-06 




0.1444 


3.0e-04 


6.8e-05 


9.8e-06 


J 


0.001205 


1.3e-04 


2.6e-05 


3.7e-06 


Zp 


0.82865 


5.2e-05 


1.6e-05 


1.9e-06 


GRV2 




1.2e-04 


1.8e-05 


5.0e-06 


GRV3 




1.6e-04 


2.4e-05 


6.5e-06 


\l-M m /M oat \ 




9.1e-05 


3.5e-05 


9.6e-06 


1 1 — ■An/'AjLitl 




4.6e-04 


8.4e-05 


2.4e-05 



Table 6. Results for a strange star model at the mass-shedding limit with p c — 3. For the meaning 
of the quantities, see Table 3. 

expect a failure of the spectral methods. Nevertheless, since the singularities show up only in 
higher derivatives, it is possible to achieve a slow convergence, see tables 4-6 13 . However, it is 
then necessary to modify the coordinate transformations (11711 81 such that the curves s = const 
do not possess a cusp (except for 5=1). Here we use 



p 2 = r 2 e st, ?=s(l-t)[r* + stg(m, (s,t)el 2 (68) 
13 Note that the example listed in Table 5 has previously been calculated, see Nozawa et al. ( 1998 1. 
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Fig. 3. Meridional cross-section of a homogeneous configuration with infinite central pressure, 
specified in Table 7 (with the axes scaled identically). The dashed curves indicate the boundary 
of the corresponding ergo-region. 



and 



rh . (l-t)[rl + stg(t)] 



P 2 = ^, <T 2 = 1 , (s,t)el 2 (69) 

for the interior and exterior region respectively. 

From the numerical results listed in Tables 4-6 we may speculate that the behaviour of the 
pressure at the star's surface (which is determined by the equation of state) affects the type of the 
above singularities. They seem to be weaker for smoother equations of state, when the pressure 
and some higher derivatives vanish at the equator. 

4.2.2. Stars with infinite central pressure 

Another possible endpoint of a sequence of rotating stars in General Relativity is reached when 
the pressure diverges at the star's centre. For example, the sequence of static homogeneous con- 
figurations is characterized by this limit. Here, the star is spherical (r p = r e and g = 0), and 
the corresponding gravitational fields are analytically given by the Schwarzschild solution which 
reads in our chosen coordinates (with r 2 - p 2 + £ 2 ) 
v , = 3 [1 - M/(2r e )] _ \-Mr 2 l{2rl) 
2[l+M/(2r e )] 2 + Mr 2 /(r 3 e ) 



_ v _ [l+M/(2r e )] 3 
1 +Mr 2 l{2rl) 

(72) 



Ws =P \ . TTTTTrk ( 71 ) 



inside (i. e. for r < r e ) and 

6 = TTm7(27) (73) 

We~ v = p [1 +M/(2r)] 2 (74) 

co = (75) 
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m=10 


m=16 


m=22 


Pc 




CO 








r p /r e 




0.7 








Q 




1.765 


4.4e-04 


8.0e-06 


1.0e-06 


M 




0.1804 


6.7e-04 


1.7e-04 


2.5e-05 






0.3865 


4.3e-04 


1.3e-04 


2.0e-05 


/ 




0.02984 


2.1e-03 


4.0e-04 


6.0e-05 






5.179 


2.9e-03 


5.0e-04 


7.5e-05 


|1 -M s 


JM mt \ 




1.6e-04 


9.7e-05 


4.3e-05 


\l-Jm 


1 ^outl 




4.7e-04 


1.6e-04 


6.8e-05 



Table 7. Results for a homogeneous model with infinite central pressure and r p /V e = 0.7. For the 
meaning of the quantities, see Table 1 . The virial identities are not defined for p c — oo since the 
integrals I\,l2,J\, J 2 diverge. 



outside the star (r > r e ). In the limit Mr" 1 — > 1 the central value e^ vanishes which corresponds 
to p c — > oo since the surface potential Vq = U'(r = r e ) remains finite. 

A rotating configuration with an infinite central pressure is characterized by an ergo-region 
that extends in the inside up to the coordinate origin, see Fig. 3. Hence, at this point the space-time 
violates the requirement of stationarity 14 and therefore some irregular behaviour of the gravita- 
tional potentials arises here, which, in the slow rotation limit, has been studied by Chandrasekhar 



& Miller ( 1974). Consequently, we again expect a failure of the AKM-method. But as in the case 
when the mass-shedding limit occurs, we are still able to obtain slowly converging numerical 
solutions, see Table 7. It is however necessary (i) to modify the Chebyshev representation of the 
interior gravitational potentials and (ii) to introduce a slightly different coordinate mapping of 
the interior region. In particular, this is done by writing 

e u ' = e v ° +0- l)Hu,(s,t) (76) 

a'e u ' = p 2 [a' B (t) + (s-l)H a ,(s,t)] (77) 

W {iM) e- u ' = p e- y ° [l + W B (t) + (s- l)H WM (s, t)] (78) 

and using again the transformation J65I . The reformulation of the interior Chebyshev expansions 
is motivated by the Schwarzschild solution given above. We learn from here that the vanishing 
of e u ' c coincides with that of the central value of W/p. Moreover, we note that the combination 

P~ 2 g w = (p-'We-"') 2 - {p- l a't u ' f (79) 



14 A locally stationary, axisymmetric spacetime requires the existence of a timelike linear combination 
of the two Killing vectors corresponding to stationarity and axisymmetry. The latter one vanishes on the 
symmetry axis. 
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9.4e-08 


6.8e-10 
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8.371248e-04 


4.0e-05 
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1.027320e-01 


5.6e-06 


2.2e-07 


2.5e-09 


./ 


3.703716e-06 


1.7e-04 


2.7e-06 


3.7e-08 


z P 


1. 59297 le-02 


3.3e-05 


1.3e-06 


1.6e-08 


GRV2 




4.2e-06 


4.7e-08 


6.3e-10 


GRV3 




6.0e-06 


6.2e-08 


8.5e-10 


|1 - M in /M out | 




2.2e-06 


5.4e-08 


4.2e-10 


1 1 — ^in/ ^outl 




1.0e-04 


3.8e-07 


1.8e-09 



Table 8. Results for a homogeneous model with p c - 0.002 and r p /r e - 0.2. For the meaning of 
the quantities, see Table 1 . 




Fig. 4. Meridional cross-section of a highly distorted homogeneous configuration at the mass 
shedding limit, specified in Table 9 (with the axes scaled identically). 



remains positive (and finite) at the origin when & u ° — > 0. So the above reformulation ensures 
particular dependencies of the metric functions at the origin when t u ' c — > 0. 

The use of the transformation \65\ allows one to lay the coordinate mesh more densely in the 
vicinity of the origin. This helps to take the singular behaviour in higher derivatives of the func- 
tions Hjji ,H a > and ff^jnt into account and thus provides a better convergence. For the example 
given in Table 7 and Fig. 3, we used c s = 0.65. In the approximation scheme we prescribed the 
parameters (e u ' c ,r p /r e ) and finally pushed e^ to zero. 
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Lle-02 


1.3e-03 


1.3e-04 
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0.9883 


1.5e-03 
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1.9e-06 
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0.03900 


6.9e-03 


9.8e-04 


5.5e-05 




0.4005 


2.5e-03 


6.9e-05 


2.3e-05 
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0.002717 


1.2e-02 


1.7e-03 


9.5e-05 


z P 


0.2388 


4.8e-03 


7.8e-04 


4.4e-05 


GtfV2 




6.2e-05 


2.3e-05 


5.9e-06 


GflV3 




1.3e-04 


3.0e-05 


7.3e-06 


\l-MJM mt \ 




1.5e-04 


2.1e-05 


4.6e-06 


1 1 — -An/«A>utl 




2.0e-03 


1.3e-05 


8.8e-06 



Table 9. Results for a highly distorted homogeneous configuration at the mass shedding limit 
with p c = 0.003. For the meaning of the quantities, see Table 1. 

4.2.3. Highly flattened stars 

As with the situations above, the coordinate transformations dl7ll8t need to be modified if one 
wants to calculate models of strongly distorted stars (such as the examples given in table 2 
of Ansorg et al. 2002). When using J17I1 8i . then each curve s = const represents an image 
which is similar to the star's boundary. This leads for distorted stars to a non-uniform partition 
of the domains by the coordinate net of (s, f)-variables. Moreover, the oblateness of the config- 
urations suggests adapting the coordinates s and t for the exterior domain to resemble oblate 
spheroidal coordinates. A possible mapping which takes these considerations into account is 
given by [(s, t) e I 2 ]: 

P 2 = rltT(t), (80) 
e = s(l-t)[r 2 p + tg(t)] (81) 

for the interior region and 

p 2 = t[rl - r\ + f(s)][l - s + st(*)], (82) 
f = (l-t)[f-(s) + tg(t)] (83) 

exterior to the star with 

1 - c t (l - t) s 
For the examples (b) and (c) in Ansorg et al. ( 2002 1 we took c s = 0.07 and c, = -0.2 whereas for 



the example (a) in Ansorg et al. ( 2002 1 as well as for that presented in Table 8 the values c s = 0.2 
and c t = were chosen. Note that the prescribed parameter pair, (p c , r p /r e ), is the same for the 
latter two configurations. All other physical quantities are also very similar for these models, 
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including the appearance of the corresponding cross sections. Here, the high accuracy is needed 
in order to distinguish between these two nearby configurations. 

A final example (see Table 9 and Figure 4) exhibits that even in the highly flattened regime, 
configurations at the mass shedding limit can be calculated (here c s = 0.2 and c t = 0). This par- 
ticular model is close to the configuration (a) in Ansorg et al. ( 2003a I, which marks the transition 
body between spheroidal and toroidal configurations at the mass shedding limit. A more detailed 
investigation of highly flattened homogeneous bodies in General Relativity will be published 
elsewhere. 
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